package fem2;

import inf.math.LinAlg;
import math2.IsoparametricEdgeGeometry;
import fem2.enu.MeshPartType;

/**
 * Edge is defined as a straight line connect equal-space point which specify by
 * n1 n3 n4 ... n2; n1: start point n2: end point. The reason for this is to
 * make consistent with LagrangeFB definition.
 * 
 * @author hbui
 * 
 */
public class Edge extends MeshPart {

	private Node[] nodes;

	/**
	 * constructor for isoparametric edge
	 * 
	 * @param nodes
	 */
	public Edge(Node... nodes) {
		this.nodes = nodes;
		if (nodes.length == 2) {
			setType(MeshPartType.LINE2);
		} else if (nodes.length == 3) {
			setType(MeshPartType.LINE3);
		}
		setBasis(BasisManager.getInstance().getBasis(getType()));
		setGeometry(IsoparametricEdgeGeometry.getInstance());
	}

	@Override
	public Node[] getNodes() {
		return nodes;
	}

	@Override
	public double[][] gGradientsAt(double... xi) {
		// int n = getBasis().getDimension();
		// int dim = getDimension();
		// double[][] gGradient = new double[n][dim];
		// double[][] dXdxi = jacXAt(xi);
		// for (int i = 0; i < n; i++) {
		// double[] dNdxi = getBasis().getBasisFunction(i).gradientAt(xi);
		// for (int j = 0; j < dim; j++) {
		// gGradient[i][j] = dNdxi[j] / dXdxi[0][j];
		// }
		// }
		int n = getBasis().getDimension();
		int dim = getDimension();
		double[][] gGradient = new double[n][1];
		double[][] dXdxi = jacXAt(xi);
		double jac = LinAlg.norm2(dim, dXdxi[0]);
		for (int i = 0; i < n; i++) {
			double[] dNdxi = getBasis().getBasisFunction(i).gradientAt(xi);
			gGradient[i][0] = dNdxi[0] / jac;
		}
		return gGradient;
	}

	/**
	 * Evaluates the Jacobian of the solutions
	 * 
	 * <pre>
	 * | duh1/dxi duh2/dxi	.	.|
	 * </pre>
	 * 
	 * at the specified point.
	 * 
	 * @param xi
	 *            natural coordinate
	 * @return the Jacobian matrix of size 1 x n
	 */
	@Override
	public double[][] jacUAt(double... xi) {
		int d = countFields();
		double[][] dUdxi = new double[1][d];

		for (int i = 0; i < d; i++) {
			dUdxi[0][i] = getSolution(i).gradientAt(xi)[0];
		}
		return dUdxi;
	}

	@Override
	public double[][] gJacUAt(double... xi) {
		int d = countFields();
		int n = getBasis().getDimension();
		int dim = getDimension();
		double[][] jq = jacXAt(xi);
		double[][] ju = this.jacUAt(xi);
		double[][] gju = new double[dim][d];

		for (int i = 0; i < n; i++) {
			for (int j = 0; j < dim; j++) {
				gju[j][i] = ju[0][i] / jq[0][j];
			}
		}
		return gju;
	}

	public double getLength() {
		// double[] X1 = nodes[0].getPosition();
		// double[] X2 = nodes[1].getPosition();
		// return LinAlg.diffNorm2(getDimension(), X1, X2);

		// return computeVolume();

		IntegrationManager integrator = IntegrationManager.getInstance();
		double[] xi = integrator.getIntegrationPoint(getType(), 1, 0);
		double w = integrator.getIntegrationWeight(getType(), 1, 0);
		double[] jac = jacXAt(xi)[0];
		return w * LinAlg.norm2(getDimension(), jac);
	}
}
